# Replication package for 
# "Latent Territorial Threat and Democratic Regime Reversals"
# Johannes Karreth, Jaroslav Tir, and Douglas M. Gibler
# Forthcoming in the Journal of Peace Research, 2021
# jkarreth@ursinus.edu

# Please read "0_readme.html" before using this code.

# This file is `4_analyze_reversals1816.R`
# It generates the results in Table A8 (Territorial threat and democratic reversals, 1816-2016).

rm(list = ls())

####################################
### Merge all datasets
####################################

set.seed(123)
library("tidyverse")

load(file = "Source/polity.RData")
polity.dat <- polity.dat[polity.dat$ccode != 769, ]


# Fix Germany int polity.dat to not lose the pre-1990 observations
polity.dat[polity.dat$ccode == 260 & polity.dat$year >= 1955 & polity.dat$year <= 1990, ]$ccode <- 255

# Territorial threat: estimated models using DR.TTpred.R, then saved as TTfm_mcmc.RData and TTpred.csv
tt <- read.csv(file  = "Work/TTpred.csv")

# Fix Germany to not lose the pre-1990 observations
tt[tt$ccode == 260 & tt$year >= 1955 & tt$year <= 1990, ]$ccode <- 255

# Create rolling mean (previous 2 & 5 years) for TT

# moving mean for that year and previous year (e.g. 5 represents the mean of that year and the four previous years)
tt2 = tt %>%
  group_by(ccode) %>%
  arrange(ccode, year) %>%
  mutate(max_terrthreatfm_mcpp_cur2 = zoo::rollmean(x = max_terrthreatfm_mcpp, 2, align = "right", fill = NA),
         max_terrthreatfm_mc1946pp_cur2 = zoo::rollmean(x = max_terrthreatfm_mc1946pp, 2, align = "right", fill = NA),
         max_terrthreatfm_mcpp_cur5 = zoo::rollmean(x = max_terrthreatfm_mcpp, 5, align = "right", fill = NA),
         max_terrthreatfm_mc1946pp_cur5 = zoo::rollmean(x = max_terrthreatfm_mc1946pp, 5, align = "right", fill = NA))
head(tt2, 75)

# moving mean for the previous years not including the current year (e.g. 5 represents the mean of the 5 previous years)
tt2 = tt2 %>%
  mutate(max_terrthreatfm_mcpp_lag1 = lag(max_terrthreatfm_mcpp, n = 1),
         max_terrthreatfm_mc1946pp_lag1 = lag(max_terrthreatfm_mc1946pp, n = 1)) %>%
  mutate(max_terrthreatfm_mcpp_pre2 = zoo::rollapply(data = max_terrthreatfm_mcpp_lag1, 
                                     width = 2, 
                                     FUN = mean, 
                                     align = "right", 
                                     fill = NA, 
                                     na.rm = TRUE),
         max_terrthreatfm_mc1946pp_pre2 = zoo::rollapply(data = max_terrthreatfm_mc1946pp_lag1, 
                                                     width = 2, 
                                                     FUN = mean, 
                                                     align = "right", 
                                                     fill = NA, 
                                                     na.rm = TRUE),
         max_terrthreatfm_mcpp_pre5 = zoo::rollapply(data = max_terrthreatfm_mcpp_lag1, 
                                                width = 5, 
                                                FUN = mean, 
                                                align = "right", 
                                                fill = NA, 
                                                na.rm = TRUE),
         max_terrthreatfm_mc1946pp_pre5 = zoo::rollapply(data = max_terrthreatfm_mc1946pp_lag1, 
                                                     width = 5, 
                                                     FUN = mean, 
                                                     align = "right", 
                                                     fill = NA, 
                                                     na.rm = TRUE))

head(tt2, 75)

# for those obs where no 2-year rolling average could be calculated, use the current value instead
tt2$max_terrthreatfm_mcpp_pre2 <- ifelse(is.na(tt2$max_terrthreatfm_mcpp_pre2) == TRUE, tt2$max_terrthreatfm_mcpp, tt2$max_terrthreatfm_mcpp_pre2)
tt2$max_terrthreatfm_mc1946pp_pre2 <- ifelse(is.na(tt2$max_terrthreatfm_mc1946pp_pre2) == TRUE, tt2$max_terrthreatfm_mc1946pp, tt2$max_terrthreatfm_mc1946pp_pre2)
tt2$max_terrthreatfm_mcpp_pre5 <- ifelse(is.na(tt2$max_terrthreatfm_mcpp_pre5) == TRUE, tt2$max_terrthreatfm_mcpp, tt2$max_terrthreatfm_mcpp_pre5)
tt2$max_terrthreatfm_mc1946pp_pre5 <- ifelse(is.na(tt2$max_terrthreatfm_mc1946pp_pre5) == TRUE, tt2$max_terrthreatfm_mc1946pp, tt2$max_terrthreatfm_mc1946pp_pre5)

dat <- left_join(x = polity.dat, y = tt2[ , c("ccode", "year","max_terrthreatfm_mcpp", "max_terrthreatfm_mcppSD", "max_terrthreatfm_mcpp10", "max_terrthreatfm_mcpp90", "max_terrthreatfm_mcpp_cur2", "max_terrthreatfm_mcpp_cur5", "max_terrthreatfm_mcpp_lag1", "max_terrthreatfm_mcpp_pre2", "max_terrthreatfm_mcpp_pre5", "max_terrthreatfm_mc1946pp", "max_terrthreatfm_mc1946ppSD", "max_terrthreatfm_mc1946pp10", "max_terrthreatfm_mc1946pp90", "max_terrthreatfm_mc1946pp_cur2", "max_terrthreatfm_mc1946pp_cur5", "max_terrthreatfm_mc1946pp_lag1", "max_terrthreatfm_mc1946pp_pre2", "max_terrthreatfm_mc1946pp_pre5"), ], by = c("ccode" = "ccode", "year" = "year"), all.x = TRUE, all.y = FALSE)

dat <- arrange(dat, ccode, year)
dat <- group_by(dat, ccode)

# Transformations of TT estimate (based on 1816-2016 sample)

# Deciles - from the whole population (not just the reversal sample)
dat$max_terrthreatfm_mcpp_q <- cut(dat$max_terrthreatfm_mcpp, breaks = quantile(dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q <- as.numeric(dat$max_terrthreatfm_mcpp_q)
dat$max_terrthreatfm_mcpp_q_pre2 <- cut(dat$max_terrthreatfm_mcpp_pre2, breaks = quantile(dat$max_terrthreatfm_mcpp_pre2, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q_pre2 <- as.numeric(dat$max_terrthreatfm_mcpp_q_pre2)
dat$max_terrthreatfm_mcpp_q_pre5 <- cut(dat$max_terrthreatfm_mcpp_pre5, breaks = quantile(dat$max_terrthreatfm_mcpp_pre5, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q_pre5 <- as.numeric(dat$max_terrthreatfm_mcpp_q_pre5)
dat$max_terrthreatfm_mcpp.lag_q <- cut(dat$max_terrthreatfm_mcpp.lag, breaks = quantile(dat$max_terrthreatfm_mcpp.lag, probs = seq(0, 1, by = 0.1), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp.lag_q <- as.numeric(dat$max_terrthreatfm_mcpp.lag_q)
dat$max_terrthreatfm_mcpp_q5 <- cut(dat$max_terrthreatfm_mcpp, breaks = quantile(dat$max_terrthreatfm_mcpp, probs = seq(0, 1, by = 0.05), na.rm = TRUE), include.lowest = TRUE)
dat$max_terrthreatfm_mcpp_q5 <- as.numeric(dat$max_terrthreatfm_mcpp_q5)

# Log

dat$max_terrthreatfm_mcpp_ln <- log(dat$max_terrthreatfm_mcpp + 0.00025)
dat$max_terrthreatfm_mcppSD_ln <- log(dat$max_terrthreatfm_mcppSD + 0.00025)
dat$max_terrthreatfm_mcpp10_ln <- log(dat$max_terrthreatfm_mcpp10 + 0.00025)
dat$max_terrthreatfm_mcpp90_ln <- log(dat$max_terrthreatfm_mcpp90 + 0.00025)
dat$max_terrthreatfm_mcpp_pre2_ln <- log(dat$max_terrthreatfm_mcpp_pre2 + 0.00025)
dat$max_terrthreatfm_mcpp_pre5_ln <- log(dat$max_terrthreatfm_mcpp_pre5 + 0.00025)

dat <- mutate(dat, aut6tran.reg.lag = lag(aut6tran.reg, n = 1), aut7tran.reg.lag = lag(aut7tran.reg, n = 1), polity2.lag = lag(polity2, n = 1), max_terrthreatfm_mcpp_ch1 = max_terrthreatfm_mcpp - lag(max_terrthreatfm_mcpp, n = 1), max_terrthreatfm_mcpp_ch2 = max_terrthreatfm_mcpp - lag(max_terrthreatfm_mcpp, n = 2), max_terrthreatfm_mcpp.lag = lag(max_terrthreatfm_mcpp, n = 1))

dat$max_terrthreatfm_mcpp_pre2_ln_std <- as.numeric(scale(dat$max_terrthreatfm_mcpp_pre2_ln, center = TRUE, scale = TRUE))
dat$aut6tran.reg.lag_std <- as.numeric(scale(dat$aut6tran.reg.lag, center = TRUE, scale = TRUE))
dat$dem6prop.global_std <- as.numeric(scale(dat$dem6prop.global, center = TRUE, scale = TRUE))
dat$polity2.lag_std <- as.numeric(scale(dat$polity2.lag, center = TRUE, scale = TRUE))
dat$aut6tran.sum_std <- as.numeric(scale(dat$aut6tran.sum, center = TRUE, scale = TRUE))
dat$yasdem6_ln_std <- as.numeric(scale(log(dat$yasdem6), center = TRUE, scale = TRUE))


dat <- filter(dat, !is.na(reversal6sample) & reversal6sample == 1)

library("rstanarm")

main1_1816mcmc_empty <- stan_glm(aut6tran ~ max_terrthreatfm_mcpp_pre2_ln,
                         family = binomial(link = "logit"),
                         data = dat, 
                         prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
                         prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
                         cores = parallel::detectCores(),
                         iter = 5000,
                         control = list(adapt_delta = 0.99))

main1_1816mcmc <- stan_glm(aut6tran ~ max_terrthreatfm_mcpp_pre2_ln + dem6prop.global_std + polity2.lag_std + yasdem6_ln_std,
                         family = binomial(link = "logit"),
                         data = dat, 
                         prior = cauchy(location = 0, scale = 10, autoscale = TRUE),
                         prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
                         QR = TRUE,
                         cores = parallel::detectCores(),
                         iter = 5000,
                         control = list(adapt_delta = 0.99))

main2_1816mcmc_empty <- stan_glm(aut6tran ~ max_terrthreatfm_mcpp_q_pre2,
                         family = binomial(link = "logit"),
                         data = dat, 
                         prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
                         prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
                         cores = parallel::detectCores(),
                         iter = 5000,
                         control = list(adapt_delta = 0.99))

main2_1816mcmc <- stan_glm(aut6tran ~ max_terrthreatfm_mcpp_q_pre2 + dem6prop.global_std + polity2.lag_std + yasdem6_ln_std,
                         family = binomial(link = "logit"),
                         data = dat, 
                         prior = cauchy(location = 0, scale = 2.5, autoscale = TRUE),
                         prior_intercept = cauchy(location = 0, scale = 10, autoscale = TRUE),
                         QR = TRUE,
                         cores = parallel::detectCores(),
                         iter = 5000,
                         control = list(adapt_delta = 0.99))

library("BayesPostEst")

mcmcReg(list(main1_1816mcmc_empty, main1_1816mcmc, main2_1816mcmc_empty, main2_1816mcmc),
        custom.model.names = c("TT logged", "TT logged", "TT deciles", "TT deciles"),
        coefnames = list(c("Intercept", "Territorial threat (logged)"),
                         c("Intercept", "Territorial threat (logged)", "Perc. democratic (global)", "Polity (t-1)", "Years as democracy (logged)"),
                         c("Intercept", "Territorial threat (deciles)"),
                         c("Intercept", "Territorial threat (deciles)", "Perc. democratic (global)", "Polity (t-1)", "Years as democracy (logged)")),
        ci = 0.90,
        hpdi = FALSE,
        caption = "Posterior estimates: Territorial threat and democratic reversals, 1816-2016.",
        caption.above = TRUE,
        format = "tex",
        file = "Tables/tab_main_1816")

# xtable version

main1_1816mcmc_empty_tab <- rbind(mcmcTab(main1_1816mcmc_empty, Pr = TRUE), data.frame(matrix(NA, nrow = 3, ncol = 6, dimnames = list(NULL, colnames(mcmcTab(main1_1816mcmc_empty, Pr = TRUE))))))
main1_1816mcmc_tab <- mcmcTab(main1_1816mcmc, Pr = TRUE)
main2_1816mcmc_empty_tab <- rbind(mcmcTab(main2_1816mcmc_empty, Pr = TRUE), data.frame(matrix(NA, nrow = 3, ncol = 6, dimnames = list(NULL, colnames(mcmcTab(main2_1816mcmc_empty, Pr = TRUE))))))
main2_1816mcmc_tab <- mcmcTab(main2_1816mcmc, Pr = TRUE)

rownames(main1_1816mcmc_empty_tab) <- rownames(main1_1816mcmc_tab)
rownames(main2_1816mcmc_empty_tab) <- rownames(main1_1816mcmc_tab)
rownames(main2_1816mcmc_tab) <- rownames(main1_1816mcmc_tab)

main_1816ft <- cbind(main1_1816mcmc_empty_tab, main1_1816mcmc_tab, main2_1816mcmc_empty_tab, main2_1816mcmc_tab)

main_1816ft <- data.frame(main_1816ft)
main_1816ft <- add_row(.data = main_1816ft, Median = as.character(nrow(model.matrix(main1_1816mcmc_empty))), Median.1 = as.character(nrow(model.matrix(main1_1816mcmc))), Median.2 = as.character(nrow(model.matrix(main2_1816mcmc_empty))), Median.3 = as.character(nrow(model.matrix(main2_1816mcmc))))
main_1816ft <- add_row(.data = main_1816ft, Median = "TT logged", Median.1 = "TT logged", Median.2 = "TT deciles", Median.3 = "TT deciles", .before = 1)

main_1816ft$Variable <- c("", "Intercept", "Territorial threat", "Perc. democratic (global)", "Polity (t-1)", "Years as democracy (logged)", "N (country-years)")

main_1816ft <- dplyr::select(main_1816ft, -Lower, -Upper, -Lower.1, -Upper.1, -Lower.2, -Upper.2, -Lower.3, -Upper.3)
main_1816ft$Variable.1 <- NA
main_1816ft$Variable.2 <- NA
main_1816ft$Variable.3 <- NA

names(main_1816ft) <- c("", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)", "", "Median", "Std. dev.", "Pr(Estimate)")

library("xtable")

main_1816xt <- xtable(main_1816ft, caption = "Posterior estimates: Territorial threat and democratic reversals, 1816-2016.",
                     label = "tab_main_1816", digits = 2,
                     align = "llccccccccccccccc")

print(main_1816xt, file = "Tables/tab_main1816.tex",
      include.colnames = TRUE, include.rownames = FALSE,
      booktabs = TRUE, caption.placement = "top", size = "tiny", digits = 2)

#################################################
# End of script
#################################################

# > sessionInfo()
# R version 4.0.3 (2020-10-10)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Big Sur 10.16
# 
# Matrix products: default
# LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] dclone_2.3-0       Matrix_1.2-18      coda_0.19-4        knitr_1.30         pastecs_1.3.21     cobalt_4.2.4       Matching_4.9-7    
# [8] separationplot_1.3 foreign_0.8-80     MASS_7.3-53        Hmisc_4.4-2        Formula_1.2-4      survival_3.2-7     lattice_0.20-41   
# [15] RColorBrewer_1.1-2 cowplot_1.1.0      BayesPostEst_0.2.1 rstanarm_2.21.1    Rcpp_1.0.6         xtable_1.8-4       countrycode_1.2.0 
# [22] reshape_0.8.8      forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4        readr_1.4.0        tidyr_1.1.2       
# [29] tibble_3.0.4       ggplot2_3.3.2      tidyverse_1.3.0   
# 
# loaded via a namespace (and not attached):
#   [1] utf8_1.1.4           tidyselect_1.1.0     lme4_1.1-26          htmlwidgets_1.5.3    grid_4.0.3           miscTools_0.6-26    
# [7] jtools_2.1.1         munsell_0.5.0        codetools_0.2-16     statmod_1.4.35       DT_0.16              miniUI_0.1.1.1      
# [13] withr_2.3.0          colorspace_2.0-0     highr_0.8            rstudioapi_0.13      ROCR_1.0-11          stats4_4.0.3        
# [19] gbRd_0.4-11          bayesplot_1.7.2      Rdpack_2.1           labeling_0.4.2       rstan_2.21.2         mnormt_2.0.2        
# [25] clarkeTest_0.1.0     farver_2.0.3         vctrs_0.3.6          generics_0.1.0       xfun_0.20            R6_2.5.0            
# [31] markdown_1.1         VGAM_1.1-4           bitops_1.0-6         assertthat_0.2.1     promises_1.1.1       scales_1.1.1        
# [37] nnet_7.3-14          texreg_1.37.5        gtable_0.3.0         processx_3.4.5       sandwich_3.0-0       rlang_0.4.10        
# [43] splines_4.0.3        checkmate_2.0.0      rjags_4-10           broom_0.7.3          inline_0.3.17        yaml_2.2.1          
# [49] reshape2_1.4.4       abind_1.4-5          modelr_0.1.8         unmarked_1.0.1       threejs_0.3.3        crosstalk_1.1.0.1   
# [55] backports_1.2.1      httpuv_1.5.4         rsconnect_0.8.16     DAMisc_1.6.2         tools_4.0.3          psych_2.0.12        
# [61] runjags_2.0.4-6      optiscale_1.2        ellipsis_0.3.1       raster_3.4-5         R2WinBUGS_2.1-21     ggridges_0.5.2      
# [67] plyr_1.8.6           base64enc_0.1-3      ps_1.5.0             prettyunits_1.1.1    rpart_4.1-15         zoo_1.8-8           
# [73] cluster_2.1.0        haven_2.3.1          fs_1.5.0             survey_4.0           magrittr_2.0.1       data.table_1.13.6   
# [79] openxlsx_4.2.3       colourpicker_1.1.0   lmtest_0.9-38        reprex_0.3.0         effects_4.2-0        tmvnsim_1.0-2       
# [85] matrixStats_0.57.0   hms_0.5.3            shinyjs_2.0.0        mime_0.9             evaluate_0.14        shinystan_2.5.0     
# [91] rio_0.5.16           jpeg_0.1-8.1         readxl_1.3.1         gridExtra_2.3        rstantools_2.1.1     compiler_4.0.3      
# [97] bdsmatrix_1.3-4      V8_3.4.0             crayon_1.3.4         minqa_1.2.4          StanHeaders_2.21.0-7 htmltools_0.5.0     
# [103] later_1.1.0.1        RcppParallel_5.0.2   lubridate_1.7.9.2    DBI_1.1.0            ggmcmc_1.5.0         dbplyr_2.0.0        
# [109] boot_1.3-25          AICcmodavg_2.3-1     car_3.0-10           cli_2.2.0            mitools_2.4          rbibutils_2.0       
# [115] insight_0.11.1       igraph_1.2.6         pkgconfig_2.0.3      sp_1.4-4             R2jags_0.6-1         xml2_1.3.2          
# [121] dygraphs_1.1.1.6     plm_2.2-5            rvest_0.3.6          snakecase_0.11.0     callr_3.5.1          digest_0.6.27       
# [127] janitor_2.0.1        rmarkdown_2.6        cellranger_1.1.0     htmlTable_2.1.0      maxLik_1.4-6         curl_4.3            
# [133] shiny_1.5.0          gtools_3.8.2         nloptr_1.2.2.2       lifecycle_0.2.0      nlme_3.1-149         jsonlite_1.7.2      
# [139] carData_3.0-4        fansi_0.4.1          pillar_1.4.7         GGally_2.0.0         loo_2.4.1            fastmap_1.0.1       
# [145] httr_1.4.2           pkgbuild_1.2.0       glue_1.4.2           xts_0.12.1           zip_2.1.1            png_0.1-7           
# [151] shinythemes_1.1.2    pander_0.6.3         stringi_1.5.3        caTools_1.18.0       latticeExtra_0.6-29 
